#coding:utf8
#https://www.cnblogs.com/qw12/p/6079244.html
#矩阵的正交三角分解（QR分解）
from fractions import Fraction

import numpy as np

def gram_schmidt(A):
    """Gram-schmidt正交化"""
    Q=np.zeros_like(A)
    cnt = 0
    for a in A.T:
        u = np.copy(a)
        for i in range(0, cnt):
            u -= np.dot(np.dot(Q[:, i].T, a), Q[:, i]) # 减去待求向量在以求向量上的投影
        e = u / np.linalg.norm(u)  # 归一化
        Q[:, cnt] = e
        cnt += 1
    R = np.dot(Q.T, A)
    return (Q, R)

def givens_rotation(A):
    """Givens变换"""
    (r, c) = np.shape(A)
    Q = np.identity(r)
    R = np.copy(A)
    (rows, cols) = np.tril_indices(r, -1, c)
    for (row, col) in zip(rows, cols):
        if R[row, col] != 0:  # R[row, col]=0则c=1,s=0,R、Q不变
            r_ = np.hypot(R[col, col], R[row, col])  # d
            c = R[col, col]/r_
            s = -R[row, col]/r_
            G = np.identity(r)
            G[[col, row], [col, row]] = c
            G[row, col] = s
            G[col, row] = -s
            R = np.dot(G, R)  # R=G(n-1,n)*...*G(2n)*...*G(23,1n)*...*G(12)*A
            Q = np.dot(Q, G.T)  # Q=G(n-1,n).T*...*G(2n).T*...*G(23,1n).T*...*G(12).T
    return (Q, R)
#求向量x和单位向量e的Householder矩阵
# x=np.array([-2,2,1]).reshape(-1,1)#创建列向量
# e=np.array([1,0,0]).reshape(-1,1)#创建列向量
#a=||x||
#w=x-ae/||x-ae||
#H=I-2wwH
# H=householder_reflection_vector(x,e)
def householder_reflection_vector(x,e):
    a = np.linalg.norm(x)
    print("a:",a)
    w=np.array(x-a*e)
    # print("w:",w)
    w_one=w/np.linalg.norm(w)
    # print("w_one:",w_one)
    I=np.eye(len(x))#注意I为二维单位矩阵
    # print("I:",I)
    H=I-2*np.dot(w_one,w_one.T.conjugate())
    # print("H:",H)
    return H

def householder_reflection(A):
    """Householder变换"""
    (r, c) = np.shape(A)
    Q = np.identity(r)
    R = np.copy(A)
    for cnt in range(r - 1):
        x = R[cnt:, cnt]
        e = np.zeros_like(x)
        e[0] = np.linalg.norm(x)
        u = x - e
        v = u / np.linalg.norm(u)
        Q_cnt = np.identity(r)
        Q_cnt[cnt:, cnt:] -= 2.0 * np.outer(v, v)
        R = np.dot(Q_cnt, R)  # R=H(n-1)*...*H(2)*H(1)*A
        Q = np.dot(Q, Q_cnt)  # Q=H(n-1)*...*H(2)*H(1)  H为自逆矩阵
    return (Q, R)

# np.set_printoptions(precision=4, suppress=True)
np.set_printoptions(formatter={'all':lambda x: str(Fraction(x).limit_denominator())})#设置np.array用分数表示
A = np.array([[0,4,1],[1,1,1],[0,3,2]],dtype=float)
# A=np.array([[1,0,1,0],[0,1,1,-1]],dtype=float)
# A = np.array([-2,2,1],dtype=float)

# (Q, R) = gram_schmidt(A)
# print(Q)
# print(R)
# print (np.dot(Q,R))
#
# (Q, R) = givens_rotation(A)
# print(Q)
# print(R)
# print( np.dot(Q,R))

# (Q, R) = householder_reflection(A)
# print("Q:",Q)
# print("R:",R)
# print ("A=QR",np.dot(Q,R))

x=np.array([4,4,2]).reshape(-1,1)#创建列向量
e=np.array([1,0,0]).reshape(-1,1)#创建列向量
# print(x,"\n",e)
H=householder_reflection_vector(x,e)
print("向量x和单位向量e的H:",H)

